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ABSTRACT 


This paper is a survey of the current knowledge available for design- 
ing and predicting the effectiveness of controllers for dyn2imic systems 
which can be modeled by ordinary differential equations. 

A short discussion of feedback control is followed by a description 
of deterministic controller design and the concept of system "state." 

Need for more realistic disturbance models led to use of stochastic pro - 
cess concepts , in particular the Gauss-Markov process. Kolmogorov and 
Wiener showed how to estimate the state and the mean-square deviation of 
the state of a stationary stochastic process from noisy measurements . 
Kalman. and Bucy extended this theory to nonstationary stochastic processes 
and gave a simpler technique for designing the estimator (or filter) that 
is well suited for use with digital computers. Feeding back this esti- 
mated state with a deterministic "no-memory" controller constitutes 
a "compensator ," a concept familiar in classical controller design. 
The system, controlled by this compensator, with random forcing functions, 
random errors in the measurements, and random initial conditions, is 
itself a Gauss-Markov random process; hence the mean-square behavior of 
the controlled system is readily predicted. As an example, a compensator 
is designed for a helicopter to maintain it in hover in a gusty wind over 


a point on the ground. 
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I. INTRODUCTION 


Automatic control is used in connection with goal-oriented, man-made 
systems . There are several, usually overlapping, reasons for using auto- 
matic control: (a) to relieve human operators from tedious, repetitive 

tasks (e.g. , dial telephones, traffic signals, automatic elevators, auto- 
pilots), (b) to speed up and/or lower costs of production processes (e.g., 
chemical process control, automatic cutting and drilling machines), and 
(c) to control rapidly-changing or complicated systems accurately and 
safely (e.g., guidance and control of spacecraft, automatic landing of 
aircraft ) . 

There are two distinctly different types of automatic control: 
programmed (or open-loop) control and feedback (or closed-loop) control 
(see Figure 1). The difference is that the output affects the signals to 
the actuator when using feedback control but does not affect them when 
using programmed control. Typical programmed control systems are those 
used with traffic signals or automatic washing machines where a sequence 
of events is scheduled by a clock-timer. Examples of simple feedback con- 
trol systems are the temperature controllers used in buildings, ovens, and 
hot water heaters, where the temperature error (the difference between 
sensed and desired temperature) is used with some controller logic, to 
switch on and off the heating element. 

The purpose of feedback control is to reduce to acceptable levels the 
effects of unpredictable disturbances and uncertain parameters, using 
acceptable levels of control action. Thus, it is important for the 
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designer of a feedback control system to have reasonably accurate models 
of (a) the system to be controlled and (b) the expected disturbances. 
Until recently, engineers had better tools for modeling the system than 
they had for modeling the expected disturbances; they used only rather 
simple descriptions such as constant values and impulse functions (or 
integrals of the latter — step functions, ramp functions, etc.). While 
such descriptions are often adequate to yield satisfactory designs when 
only one disturbance source acts on the system, more realistic descrip- 
tions are needed when several disturbance sources act on the system 
simultaneously. 

Many systems of practical interest can be modeled by linear ordinary 
differential equations. A convenient form for such a model is a set of 
coupled first-order ordinary equations: 


X = F(t)x + G(t)u + r(t)w ; x(t ) = x , 

o o 

(1.1) 

y = M(t)x , 1 = Nx(t^) , 

(1.2) 

z = H(t)x + V, 

(1.3) 


where 


x(t) 
u(t ) 
w(t) 

y(t) 

z(t) 

v(t) 

F(t) 

G(t) 


n- vector of state variables, 

m- vector of control variables, 

r- vector of process disturbance variables, 

s-vector of terminal state parameters (s < n), 

p-vector of output variables (p < n), 

q- vector of measurement variables, 

q- vector of measurement disturbance variables, 

nxn open-loop dynamics matrix, 

nxm control distribution matrix. 
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r(t) = nxr disturbance distribution matrix, 

M(t) = pxn output distribution matrix, 

N = sxn terminal state distribution matrix, and 
H(t) = qxn measurement distribution matrix. 

In most problems, the vectors are deviations from nominal values and the 
matrix elements are partial derivatives evaluated on the nominal "path. " 

The control problem is to design a dynamic system (a "compensator") 
with measurements, z , as input and controls, u , as output, so that 
y(t) and | of the controlled system are kept acceptably close to zero 
in the presence of anticipated disturbances w,v, and x^ , while using 
acceptable amounts of control, u . 

During the last ten years it has become clear that excellent compen- 
sators can be designed in two completely separate steps: (a) a state - 

feedback controller , u = -Cx , is designed to minimize an appropriate 
quadratic error criterion, assuming the availability of perfect measure- 
ments of the system state variables and no disturbances, i.e. , z = x , 
v = 0 , w = 0 ; (b) a measurement-feedback estimator (or filter ) is de- 
signed to yield maximum likelihood estimates of the state variables of 
the system, x , from the available sensor measurements, z , taking into 
accoimt w,v,x^, and assuming perfect knowledge of the control variables, 
u(t) . The state- feedback controller is then used with the state esti- 
mates from the filter, i.e., u = -C(t)x . The resulting control system 
forms a "dynamic compensator" (see Figure 2) which minimizes the quadratic 
error criterion on the average [Joseph and Tou (1961), Gunckel and 
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Franklin (1963), and Potter (1964)].* 

In Section II we discuss state feedback controller design, and in 
Section IV we discuss the design of state estimators (measurement- 
feedback filters). Section III is a review of the concept of a Gauss- 
Markov random process which forms the basis of the filter design. Sec- 
tion V discusses the behavior of systems controlled by such filter/state- 
feedback compensators. 

II. DESIGN OF STATE- FEEDBACK CONTROLLERS 

For w = V = 0 , and z = x , a rather general design problem for 
the system (1.1)-(1.2) is to try to find an acceptable u(t) program 
in some time interval, < t < t^ , so that the output variables, y(t) 
and I stay acceptably close to zero for a given . 

A useful way to judge one design relative to another is to specify 
a quadratic criterion of the form: 

J ^ ^(y^Ay + u^Bujdt (2.1) 

t 

o 

where the positive definite weighting matrices 0 , A(t) , and B(t) 

T 

are chosen to express relative preferences of the designer, and [ ] 
indicates matrix transpose of [ ], i.e. , interchange rows and columns. 

Using the calculus of variations (or dynamic programming), it is 
straightforward to show that the state-feedback controller that minimizes 

* . 

Joseph and Tou, as well as Gunckel and Franklin, showed this for Gauss- 
Markov sequences (discrete-step processes); Potter extended the proof 
to continuous Gauss-Markov processes. 
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( 2 . 2 ) 


(2.3) 


and 

S = -SF - F^S + C^BC - m\m ; S(t^) = N^)n . (2.4) 

The nonlinear equation, (2.4), for the symmetric mXn matrix, S(t) , 

T 

is called a backward matrix Ricatti equation because of the term C BC = 
-1 T 

SGB G S which is quadratic in the elements of S , and because it must 


be integrated backward from the final time, t^ , to the initial time, 

t , to generate C(t) . The significance of S is that 
o 


J . = — X (t )S(t )x(t ) . 

min 2 o o o 


(2.5) 


Examples of such time-varying terminal control problems and their solu- 
tions are given elsewhere [e.g. , Bryson and Ho (1969), Section 14,6]. 

Obviously it is not always possible to control all of the outputs 
with a given set of controls. A "controllability" criterion was given by 
Kalman, Ho, and Narendra (1963). 

Of special Interest is the class of problems with stationary 
solutions . For such problems, the matrices F,G,M,A,B must be constant, 
and t^- t must be large compared to characteristic times of the con- 
trolled system. S may then have a constant asymptotic value, which 
means that the feedback gain matrix tends to a constant value: 


♦ 

See, e.g., Bryson and Ho (1969), Section 5.2 If x(t ) is known ahead 
of time, the optimum u(t) can be pre-calculated and programmed as an 
open-loop input; this is a peculiar feature of deterministic control 
design — it predicts no difference between open-loop and closed-loop 
control. 
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t 00 . 


( 2 . 6 ) 


S(t) - Soo 
-1 T 

= B G Soo 


(2.7) 


A controller of the form (2.2) with constant gains, , is called a 

regulator . In principle, the regulator gains could be found by 

solving sets of quadratic equations [Eqs. (2.4) with S = 0 ] to find 
the elements of the steady-state matrix, . However, this turns out to 

be an impractical procedure for n > 3 . Until recently, the procedure 
most often used was to integrate the matrix Ricatti equation, (2.4), back- 
ward on a computer until S(t) reached a steady-state [see, e.g. , Kalman 
and Englar (1966)]. This turns out to be expensive and involves several 
numerical difficulties and hazards. MacFarlane (1963) and Potter (1966) 
showed another ingenious way to find using eigenvector decomposition . 

It does not seem to be widely known or used, perhaps because it requires 
an accurate computer program to find the eigenvalues and eigenvectors of 
a large (2itK2n) matrix. However, Francis (1962) developed a remarkably 
fast and accurate algorithm for eigenvalue-eigenvector determination (the 
"QR algorithm"), and Hall (1971) applied it to the MacFarlane-Potter method 
for finding Sc^ . 

The eigenvector decomposition method is based on finding the eigen- 
vectors of the Euler-Lagrange equations for the problem of minimizing 
( 2 . 1 ): 


X 


r -1 T "1 



F , -GB G 


X 

T T 



-M AM , - F 


X 

_ 


_ _ 


( 2 . 8 ) 


The eigenvalues of this system are symmetric about the imaginary axis in 
the complex plane, (if s^^ is an eigenvalue, -s^ is also an eigenvalue.) 
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Hence, the characteristic equation of (2.8) is called the symmetric- root 


characteristic equation [see Letov (1960), Chang (1961), and Rynaski and 
Whitbeck (1966)]. If we let T be the 2nXn matrix consisting of the 
n eigenvectors of (2.8) corresponding to the eigenvalues with negative 
real parts, then Potter (1966) has shown that: 

= A (X ) ^ where T = 

Furthermore, the eigenvalues of the controlled system (i.e., the eigen- 
values of F - GC) are the eigenvalues of (2.8) that have negative real 
parts, and the eigenvectors of the controlled system are the columns of X 
A revealing alternative form for (2.8), in terms of the Laplace 
transform of the control variables, u(s) , was suggested by Rynaski and 
Whitbeck (1966): 

[b + Y^(-s)AY(s)J u(s) = 0 (2.10) 

where Y(s) = M[sl - F] is the matrix of transfer functions from the 

control variables to the outputs, i.e., 

y(s) = Y(s)u(s) (2.11) 


X 

A 


(2.9) 


and s is the Laplace transform variable. The corresponding equation 
for the system using feedback control gains, C , may be written 

.-1 


I + C(sl - F)”^G u(s) = 0 . 


( 2 . 12 ) 


It follows from (2.10) and (2.12) that 


B + Y^(- s)AY(s) = I + C(- si - F)“^G 


B 


I + C(sl - F) G 


(2.13) 

which is an n^m matrix equation (where u is an m-vector). For s = ilO 
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where u) = frequency, (2.13) is the well-known "spectral factorization" 
problem encountered in solving certain Wiener-Hopf integral equations. 

The problem here is to find the mxn matrix C that satisfies (2.13), 
i.e. , to "factor" the left side of (2.13) into two mxn matrices, one 
with eigenvalues in the left half of the complex s plane and the other 
with exactly the negative of those eigenvalues which are then in the right 
half plane. 

This factorization can be done very quickly and accurately by eigen- 
vector decomposition, using (2.8) and a digital computer [see, e.g. , 

Bryson and Hall (1971)]. It can also be done as follows: (a) find the 
n eigenvalues of (2.8) that have negative real parts and the correspond- 
ing n eigenvectors of ( 2 . 10 ) (each with only m components); then (b), require 
that the eigenvalues and eigenvectors of (2.12) be the same as those in 
(a). If m = 1 , all that is needed is to equate coefficients of like 
powers of s of the characteristic equation of (2.12) and the closed-loop 
characteristic equation formed from the n eigenvalues of (2.8) that 
have negative real parts; this gives n linear equations for the n com- 
ponents of C. If m > 1 , the requirement that the eigenvectors be 
parallel adds n(m - 1) equations, giving a total of nm linear equations 
for determining the nm components of C. 

As an example , we design a regulator for a fourth-order model of a 

helicopter in hover. Figure 3 shows the nomenclature and reference axes. 

The four state variables are 6 = pitch angle of the fuselage, q = pitch 
♦ 

For m > 1, the characteristic equation of (2.10) is of order 2mn instead 
of 2n, and hence contains extraneous roots (pointed out to author by 
J. V. Breakwell). 
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rate (q = 0) , x = horizontal distance of mass center from the hover point, 
and u = horizontal velocity (u = x) . The control variable is 6 = tilt 
angle of the rotor-thrust vector with respect to the fuselage. The pri- 
mary disturbance is horizontal wind, w , which we shall consider in a 
later section. The vertical motions (controlled by increasing or decreas- 
ing rotor thrust), the lateral motions (controlled by tilting the rotor 
thrust laterally), and the yawing motions (controlled by the tail rotor) 
are all nearly uncoupled from the longitudinal motions considered here 
and may be treated as separate regulator design problems. The longitudinal 


equations of motion in state variable form may be written as: 



( 2 . 14 ) 


A reasonable quadratic performance index is: 



and A^(s) =0 is the open-loop characteristic equation. Thus (2.16) 
becomes: 


A (s)A (- s) + 
o o 


B 




s + n 




+ n 


0 . 


(2.17) 


Coefficients for a typical transport helicopter are given in Table 1, 
along with the open-loop eigenvalues and eigenvectors. Note the uncon- 
trolled helicopter is slightly unstable in hover, which is fairly typical. 
The weighting factor ^x/B may be selected as follows; rewrite (2.15) as 


J 



(2.15a) 


6 0 

5 is the maximum amount of rotor tilt we would like to com- 
mit to position control, and we would be willing to use this amount when 
the position error is 2 meters. It follows then that we should choose 

6 = 5.0/57.3 rad and x =2.0 meters. Hence 

o 




~ -2 

= .0020 m 


(2.18) 


o 

It is often helpful to the designer to see how the eigenvalues change 
as the weighting factors change. The concept of a root locus plot was 
introduced by Evans (1950); Chang (1961) applied it to the root square 
characteristic equation (2.10) and l^naski and Whitbeck (1965) extended 
the application to multi-input, multi-output systems. Figure 4 shows a 
root locus of (2.17) with ^x/B as the parameter. The open-loop poles and 
their reflections across the imaginary axis are the "poles” (shown as x's 
of the root locus; the "zeros” (shown as O’s ) are at + (.25 + 2.50j) sec 


) 

-1 
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For the weighting factor given in (2,18) the closed-loop eigenvalues 
and eigenvectors, as well as the corresponding feedback gains, are given 
in Table 1. 

Even for = 0 , (B 0) the 0,q,u states of the helicopter 
are stabilized by this technique; the two unstable complex roots of the 
open-loop system are simply "reflected" across the imaginary axis, while 
the stable real root and the zero root remain unchanged. The corresponding 
gains are: 

C. = .0535 , C„ = .0785 sec. , C = C = 0 . (2.19) 

6 q ’ u X 

III. MODELING RANDOM DISTURBANCES AS STOCHASTIC PROCESSES 

How can "unpredictable disturbances" be modeled realistically? One 
very useful way is to describe them as stochastic processes , a concept 
which was largely developed in this century by Einstein ("Brownian motion," 
1905), Langevin (1908), Markov (1912), Planck (1917), Wiener (1930), 
Uhlenbeck and Omstein (1930), Kolmogorov (1931), Khintchine (1938), 

Wold (1938), Chandrasekar (1943), Rice (1944), Wang and Uhlenbeck (1945), 
Phillips (1947), Laning and Battin (1958), Newton, Gould, and Kaiser 
(1957), and Kalman and Bucy (1961), to mention only a few. 

A special class of stochastic processes, the Gauss-Markov process 
(CMP), has emerged as a particularly simple and useful one in describing 
random disturbances. Many natural and man-made dynamic processes can be 
modeled quite well by OWP*s. A GMP may be described by giving its mean 
value vector, x(t) , and its correlation matrix C(t,T), where 
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(3.1) 


x(t) = E[x(t)] 

c(t,r) = Ej[x(t) - x(t)] [x(T) - x(r)]'^j , (3.2) 

E[ ] means "expected value of," and x is a column vector with n 
components. The operation on the right side of (3.2) is an "outer product" 
or "dyadic product"; in Cartesian tensor notation, 

C. ,(t,T) = E |x.(t) - x.(t)l Ix.(T) - x.(t)11 . 
iJ (11 1 J ' J J M 

The "expected value of [ ] " may be thought of as an "ensemble average," 

i.e., the average of [ ] over a large number of samples of the random 

process x(t) . If x^^^(t) is the i^^ sample history, then 


x(t ) 

N 

1 im 1 

(i),^x 
X (t) 

(3.3) 


i = 1 



C(t, T) 

N 

lim 1'^ 

x^^\t) - x(t)j [x^^\t) - x(T) 

(3.4) 


i = 1 




The covariance matrix , a measure, of the mean-square deviation 
of the random process from its mean value, is given by 

(3.5) 


X(t) = E|[x(t) - x(t)] [x(t) - x(t)] 


Obviously X(t) - C(t,t) . The probability density of x(t) at any 
time t is gaussian; i.e., 



where p(Od(|) is the probability that the value of x lies between 
^ and • 
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The basic GMP is the purely-random process; it is an idealized, very 


jittery process that is completely uncorrelated from one time to the next. 
Its correlation matrix has the form 

Q(t) 6(t - T) (3.7) 

where &(■) is the Dirac delta function. The purely random process may 
be thought of as the limit of a sequence of impulses, random in both magni- 
tude and time of occurrence, where the average magnitude of the impulses 
is zero (as many positive as negative impulses), but the mean square magni- 
tude is [cr(t)J . The power spectral density is then 

Q(t) = 2[o(t)]^ p(t) (3.7a) 

where 3(t) is the average number of impulses per unit time that occur at 
time t . 

A GMP can always be represented by the state vector, x(t) , of a 
linear dynamic system forced by a gausslan purely random process where the 
initial state vector is gausslan: 

X = F(t)x + r(t)w 

where 

E[w(t)] = w , E|[w(t) - w] [w(T) - w]’^| = Q(t) 6(t - T) 

(3.8) 

E[x(t^)] = x^ ,E|[x(t^) - x^] [x(t^ - x^]^| = , 

e'[ w(t) - w] [x(t ) - X = 0. 

I o O ' 

This markov representation shows clearly that the state property for 
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deterministic processes and the markov property for random processes are 

precisely analogous: (a) if knowledge of a finite-dimensional vector at 

any time, t , makes the future of a deterministic process independent of 
o 

the past, then it has the state property . For example, if w(t) in (3.8) 
is known, then knowledge of x(t^) , the initial state vector, is suffic- 
ient to predict x(t) for t > t ; (b) if knowledge of the probability 

o 

density function for a finite-dimensional vector at any time, t^ , makes 

the future of a random process Independent of the past, then it has the 

markov property . For example, if Q(t) in (3.8) is known, then knowledge 

of p[x(t )1 , the initial probability density function, is sufficient to 
o 

predict p[x(t)] for t > t^ . 

Since the gaussian probability density of a vector, x , is completely 
specified by giving its mean value x and its covariance matrix, X (see 
3.6), a GMP is specified for ^ > t^ by giving x(t^) ,X(0,F(t ) ,F(t ) , w(t ) , 
and Q(t) . In fact, it can be shown that 

x = Fx + l^;x(t)=x, (3.9) 

o o 

X = FX + XF^ + rOT^ , X(t ) = X .* (3.10) 

o o 

Note that x(t) and X(t) are completely independent of each other. It 
can also be shown (e.g. , Bryson and Ho, p. 333) that the correlation 
matrix is given by: 

! (I)(t,r) X(T) ; t > T , 

(3.11) 

X(t) $ (T,t) ; t < T , 


If F = constant matrix with eigenvalues s. , i=l,... ,n, then it 

can be shown that the ^(n + 1) eigenvalues of (3.10) are all possible 
combinations s^ + Sj , i , j = 1 , . . . , n . 
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where <|>(t,T) is the ii<n transition matrix of x = Fx , i,e., 

d){t,T) = F(t) d)(t,T) : 0(T,T) I . (3.12) 


One of the simplest nonstationary GMP's is the so-called Wiener 
process, which is the integral of stationary white noise: 

X = w 

E[x(0)] = 0 , E[x(0)]^ = 0 , (3.13) 

D[w(t)] = 0 , E[w(t) w(t)] = Q6(t -T ) , Q = constant. 



The root-mean-square (RMS) value x(t) is .yx(t ) = ./qF”, which shows 

the diffusive or random walk character of OVlP's. The Wiener process may 

be thought of as a series of random step functions of average magnitude 

zero occurring at random times but with a certain average number of steps 

2 2 

per unit time, P . It follows that Q = 2c P where c is the mean 
square magnitude of the step functions. 

Another well-known nonstationary GMP is the underdamped oscillator 
forced by white noise [Wang and Uhlenbeck (1945)]. This process could be 
interpreted, for example, as a pendulum in a random wind (see sketch 
below). The equations of motion in state variable form may be written as 
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where w(t) is proportional to wind velocity, o) = undamped natural 

frequency of the pendulum, and ^ = damping ratio. Suppose the motion 

starts with 6(0) = 6 , q(0) = 0 and w(t) = 0 , E[w(t)w(T)] = 

o 

Q 6 (t - T ), i.e. , we know the initial state exactly [Xgg(O) = X^^(O) = Xq 0 (O) = O] 
and the wind has zero mean velocity but has a constant spectral density 
proportional to Q. The mean value history is given by the familiar 
equations: 

6 = q ; 0 ( 0 ) = 6 ^ , 

q = 6 - 2 ^ 0 X 1 : q( 0 ) = 0 . (3.16) 

Tlie elements of the covariance matrix are given by 


^06 

= 2X_ 

0q 

> 

Xgg(O) = 0 , 

• 

2 

2^j) X„ + X ; 

^ 6q qq ’ 

Xg^(O) = 0 , 

^eq 

= - 

• 

2 

- 4^(C + Q ; 

qq 


X 

qq 

= - ’<8, 


For ^ =0.2 , the time histories are shown in Figure 5. Note the two 
variances and the covariance change with twice the frequency that q and 
6 change. A more revealing display of the time histories is given in 
Figure 6, where q/a> is plotted against 6 and probability ellipses are 
shown at selected points on the path. The axes of the ellipses are 
in the directions of the eigenvectors of X and the major and minor axes 
are proportional to the eigenvalues of X . Note the angular velocity q 
becomes uncertain first (elongated ellipse parallel to q axis) but the 

/ T2 

uncertainty quickly extends to 6 . Foryl - ^ cut » 2jt , 
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0 , q 


Q 

45E • 


(3.18) 


eq 


° ■ hs 


Q 


4^CD 


qq 


Thus, it is the presence of damping (^ > 0) that produces a statistical 
steady-state . The uncertainty created by the random wind is just balanced 
by the effects of damping. Note that X does not always increase with 
time; if X(0) > X(“) , then X would decrease with time toward X(oo) . 


Stationary Gauss-Markov Processes 
For a stationary GMP , 

C(t,r) = C(t - T) , (3.19) 

that is, C is a function of (t - T) only. The Fourier transform of 
the correlation matrix of a stationary GMP is called its spectral density 
matrix , Q(<jd) , i.e., 

00 

Q(cu) = ^ C(t)e'‘*^* dt (3.20) 

— 00 

which implies that 

00 

C(t) = ^ Q(co)e^^ da) (3.21) 

— 00 

The scalar, stationary, purely random GMP has a constant spectral 
density and for this reason is often called white noise, since it is made 
up of all frequencies from 0 to Its correlation is therefore 

Q6(t - T) , (3.22) 

where Q = spectral density (constant). 
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At the opposite extreme from white noise is a random bias; it is 


simply an unpredictable constant. It has a constant correlation and hence 
its spectral density is given by 

C6(cu) (3.23) 

where C/2n = correlation (constant). Its spectrum is thus a "line" at 
U) = 0 (zero frequency). In between white noise (zero correlation time), 
and a random bias (infinite correlation time), stationary GMP's have 
finite correlation times, i.e., C(t - T) — > 0 as |t - t | — > <» . For 
example, white noise forcing a first-order system produces a GMP with an 
exponential correlation (see Figure 7). As another example, the components 
of the correlation matrix for the pendulum in a random wind (when it has 
reached a statistical steady state*) are also shown in Figure 7 for ^ = 0.2 

IV. DESIGN OF STATE ESTIMATORS 

In the early 1960 *s a significant advance was made in control theory 
for random systems. Kalman and Bucy (1961), using the calculus of varia- 
tions, developed the maximum likelihood filter for estimating the state of 
a nonstationary GMP in real time from measurements of the process that con- 
tain additive white noise. This extended and clarified the pioneering work 
of Wiener (1942) and Kolmogorov (1941), who had developed the minimum mean- 
square error filter for stationary stochastic processes (not necessarily 
markovian or gaussian) using generalized harmonic analysis. 

Kalman and Bucy assumed that a mathematical model of the (MP to 
be estimated was known, (1.1), and assumed that the process noise, w(t) , 
and the measurement noise, v(t) , were purely random processes: 

it 

Note that if q < n , z(t) is, in general, not a markov process since 

z(t ) does not enable one to predict z(t) for t> t 
o o 


E[w(t)] = 0 , E|[w(t)] [wCt)]^ - Q(t) 5(t - T) 

E[v(t)] = 0 , E|[v(t)] tv(T)]r = R(t) 5(t - T) * (4.1) 

where Q(t) is a known, possibly time- varying, r X r spectral density 
matrix for the process noise and R(t) is a known, possibly time-varying, 
q X q spectral density matrix for the measurement noise. 

The Kalman-Bucy (KB) filter is then given by: 


X = Fx + Gu + K(z - 1^) 


x(t ) = x(t ) , 

o o 


(4.2) 


where 

T -1 

K = PH R 

P = FP + PF^ + rOT^ - KRK^ ; P(t^) = X(t^) 


(4.3) 


Note F,r,K,H,P,R,Q can all be time-varying quantities, and the process 
can start at any time, t^ . Gu is assumed to be a known forcing 
function. The filter (4.2) is similar to the mean value model of (3.9), 

with the addition of some known forcing functions, Gu , and some feedback 
terms proportional to (z - Hx) , the difference between the actual 
measurements and the expected value of the measurements based on the 
current estimate of the state. (The expected value of the measurement 
noise, v , is zero since a purely random process cannot be predicted. ) 

The rate of change of the estimated state, x , is thus proportional to 
the measurement prediction error. - 


^Obviously, there is no lack of generality in assuming w(t) = v(t) = 0. 
+ 

If there is uncertainty in implementing u(t), it can be included in 
w(t) 
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The errors in the estimates given by (4.2) are defined as 



By subtracting (1.1) from (4.2), a set of differential equations for the 
errors in the estimates can be found: 

X = (F - KH) X + Kv - Pw : x(t ) = x(t ) - x(t ) , (4.5) 

o o o 

From (4.5) it is clear that errors in the estimates are caused by (a) 
random errors in the measurement signals, v(t) ; (b) random forcing 
functions, w(t) , and random initial conditions, • 

The nXn symmetric matrix P(t) in (4.3) is the covariance matrix of 
the error in the estimates, i.e. , 

P(t) = Ej[x(t) - x(t)] [x(t) - x(t)]^| . (4.6) 

It may be computed from (4.3), which is a forward matrix Riccatl equation , 
T T -1 

since KRK = PH R HP is quadratic in the elements of P. Comparing 

(4.3) with (3.10) we see that they are the same except for the information 
T 

rate term (-KRK ) . Obviously P is less than or equal to X , i.e., 
with measurement information the uncertainty in the state is usually 
reduced. 

The nXq gain matrix, K(t) , is proportional to the uncertainty in 
the estimate of the state (the covariance matrix, P ) and inversely pro- 
portional to the uncertainty in the measurements (the spectral density 
matrix R ) . Thus, if the uncertainty in state is large compared to 
uncertainty in measurements, the filter weights the measurement prediction 
error heavily and vice-versa. 
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Obviously it is not always possible to estimate all the state vari- 
ables from a given set of measurements (4,1), An "observability” criterion 
was given by Kalman, Ho, and Narendra (1963), 

As an example , suppose we had a measurement of the angular velocity 
of the pendulum in a random wind, i,e, , 

z=q + v, (4,7) 


so that H = [0,1] , and assume R = constant , The KB filter is then 


e 


where 




and 



P 


q + K (z - q) ; 6(0) = e , 

C7 O 

- - 2^o)q + (z - q) ; q(0) = 0 , 


ip K 

R 0q ’ q 


- P , 

R qq ' 


2P - - 

6q R 0q 




OD^ P„„ - 2^CU P^ + P - i Pfl P ; Po (0) 
60 ^ 0q qq R 6q qq 0q 


P^ - 4t;a) P + Q - i P^ 
^ 6q ^ qq ^ R qq 


P (0) 

qq 


(4,8) 


(4.9) 



(4.10) 


The estimates 0 and q from (4.8) depend on the measurement history, 

z(t) , of the particular run. However, we can predict P(t) , the mean- 

2 

square error of the estimates, ahead of time. For ^ = 0. 2 and Q/co R = 

2.0 , the time histories of P^r, > P^ . P ^re shown in Figure 8. They 

’ 00 0q’ qq 

are similar to the time histories of X__, , X in Figure 5 but are 

00 0q qq 

less. The uncertainty created by the random wind is balanced not only by 
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damping but also by the information coming from the measurement of q . 


Again a statistical steady-state is reached for 




C^CDt 


» 2jt 



^ 0 


ee 



2 

^ + V 


p 

qq 



2 


where ^ ' 



Q 


2 

R 


(4.11) 


2 

For Q/uj R 0 (noise in the measurement » wind noise), Equations (4.11) 

2 

indicate no change from the unmeasured situation (3.18). For Q/o) R -$► oo 

(noise in the measurement « wind noise), — > 0 , P — > 0 , i.e., 

uu qq 

the errors in the estimates become negligibly small. 


Stationary Filters 


The Wiener-Kolmogorov (WK) filtering technique for stationary gaussian 
processes, while remarkably creative, was difficult to understand for 
most engineers. It involved solving an integral equation (of the Wiener- 
Hopf type) for the Fourier transform of a kernel function, W(t - T) , 
such that a scalar "signal " is estimated with minimum mean-square error, 
by t 




W(t - T) z(T) dt 


(4.12) 


where z(t) is the signal plus noise; the time correlations (hence 
spectral densities) of the signal and the noise are assumed known. Even 


* 

Their treatment did not assume that the process was Markovian. 
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after finding W(t -T) , the problem remained of implementing (4.12), 
which, being a convolution integral, is not useful for real-time compu- 
tation of the estimated signal. A dynamic system had to be found whose 
impulse-response closely matched W(t - T) . 

In engineering circles it gradually became known that almost the 
only practical way to design WK filters was to assume that both the signal 
and the noise had "rational spectra" (see, e.g., Phillips, Section 7.4, 
or Lanlng and Battin, Section 7.5), which is equivalent to assuming that 
they are components of vector GMP's. The Kalman Bucy development starts 
with this assumption and yields the dynamic system for the filter directly 
[(4.2) with all coefficient matrices constant]; finding the steady-state 
solution of (4.3) is thus equivalent to solving a Wlener-Hopf Integral 
equation for rational spectra. Unlike the WK approach, the KB approach 
handles the multiple-output, multiple-measurement case with ease. 

However, the WK filter handles measurement noise with rational 
spectra in a relatively straightforward manner, whereas the KB filter, 
as originally presented (1961), required that the measurement noise 
be white (constant spectra). Fortunately, a simple extension of the 
KB technique enables one to handle "colored" measurement noise even 
for nonstationary GMP’s [cf. Bryson and Johansen (1965)]. In this 
case the measurement appears directly as part of the estimated state 
("feedforward") and, if the measurement noise is sufficiently smooth, 
the optimum filter may even include differentiators. 

Stationary KB filter design may be done using eigenvector decompo - 
sltion. The Euler- Lagrange equations for the filter are (see Bryson 
and llo, p, 396): 
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F , -roi' X 


T -1 T 

-HR H , - F X, 


(4.13) 


T -1 
H R z 


If we let be the 2nXn matrix consisting of the n eigenvectors of 

(4.13) corresponding to the eigenvectors with positive real parts, then 


P = - X (A ) ^ where T = ^ 

a> + + + . 


(4.14) 


The eigenvalues of the error estimates, x = x - x , are the eigenvalues 
of (4.13) with negative real parts and the eigenvectors are the columns 
of [a^] ^ [see, e.g. , Bryson and Hall (1971)] . 

Again, a more revealing form for (4.13) may be given in terms of the 
Laplace transform of the measurements, z(s) , and the measurement esti- 
mates, F^(s) : 


[R + Z(s) Q Z^(-s)]R"^[ z(s) - Hx (s)] = z(s) 


(4.15) 


where Z(s) = H[sl - F] r is the matrix of transfer functions from the 
disturbance variables to the measurement variables, i.e.. 


z(s) = Z(s) w(s) . 


(4.16) 


The corresponding equation for the filter using measurement feedback 
gains, K , may be written; 


n + H(sl - F)~ K] [z- Hx(s)] = z(s) . 


(4.17) 


It follows that the spectral factorization equation is: 
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[R + Z(s)Qz’’^(- s)] = [I + H(sl - F)"^k]R[I + H(-sI - F)-^K]^ 

(4.18) 

which Is a qXq matrix equation (where z is a q-vector). The problem 
here is to find the nXq filter gain matrix K. Again, by requiring the 
eigenvalues with negative real parts of (4,13) to be the same as the 
eigenvalues of (4.17), the single measurement case is solved (q = 1) . 

For q > 1 , we require, in addition, that the eigenvectors of (4.17) be 
parallel to the eigenvectors of (4. 15) that correspond to the eigenvalues of (4.13) with 
positive real parts [see the discussion following (2.13)], 

As an example , we design a steady-state filter for the helicopter in 
hover (see the last part of Section II). We assume a position measurement 
and a pitch angle measurement, both with jittery noise having correlation 
times short compared to decay times of the closed-loop eigenvalues; thus 
they can be modeled as white noise processes. For an RMS position measure- 
ment error of 1 meter with a correlation time of .136 sec,, the spectral 
density is 

R = 2(1)^ (.136) = .272 m^ sec. (4.19) 

X 

For an RMS pitch angle measurement error of 0.5 deg. with a correlation 
time of 0.10 sec., the spectral density is 

2 

Rg = ‘^(■^- 2 ) (0-1) = 1.53 X 10 "^ (rad )^ sec. (4.20) 

On the other hand, the correlation time of wind fluctuations is comparable 
to, or larger than, the decay times of the closed- loop eigenvalues so we 
model the wind with a first-order shaping filter as 
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(4.21) 


T w + w = q , 
c w 


where q is white noise. We have taken T = 5 sec. and RMS wind 
w ^ 

magnitude as 7 msec ^ ; this implies that the spectral density of q^ 
should be 

Q = 2(7)^(5) = 490 m^ sec"^ . (4.22) 

The wind velocity, w, is now an additional (fifth) state variable. 

The state-feedback gains on d, q, u, and x are not changed by adding w 
as a state variable [in the root locus, Figure 4, a zero is added that just 


cancels the additional pole at s= -(1/T )], but an additional gain is 

c 

added for w, i.e. , 5 = -C.0 -C^-Cx-Cx-Cw.. It is straight- 

’ ’ 0 q u X w 

forward to show that 

C = b”^ S (4.23) 

w w 

where 

S = S r . 
w oo 



For the numbers used here, this gives 

C = -.00195 m ^ sec. (4.24) 

w 


The filter design, on the other hand, is significantly affected by adding 
the shaping filter. The root square characteristic equation in the form 
(4.15) is given by 

|r + Z(s)QZ^(-s)| = 0 (4.25) 

where 



Z(s) 
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z (s) ^ = L ' ""2 / S J 

x' q (s) (T s + 1) Aq(s) ’ 

w c ^ 

_ , , _ 0(s) _ 

q (s') (T s + 1) A (s) ’ 

w c o 

and given below (2.16). Thus (4.23) may be written as: 


+ Z^(s) Z.(-s) + -^ Z (s) Z (-s) = 0 . (4.25a) 

W C7 C7 K X X 

X 

A root locus of (4.25a) with Q/Rg as the parameter is shown in Figure 9 
for the numerical values of Table 1, 

1.525 X 10“® _-2 
R ~ .272 "* • 

X 

and T = .5 sec. The "zeros" of the root locus are at .190± .193j sec , 

c 

and the "poles" are the open- loop poles of Table 1, the "wind pole" at 
s = -0.2 sec and their reflections in the imaginary axis. For Q = 490 

n 

m sec“^, the eigenvalues of the filter estimate errors are given in Table 2, 
as well as at the corresponding filter gains and predicted RMS values for 
the filter estimate errors. Table 3 shows the dynamic compensator, 
comprised of the filter and the estimated state feedback to rotor tilt, 6. 

Even for Q = 0, , a stable filter is designed by this technique. 

The two unstable complex roots of the open-loop system are "reflected" 
across the imaginary axis, while the stable real roots and the zero 
root remain unchanged. The corresponding filter gains and RMS estimate 
errors are: 


r-. 00126, -.000635, 

K = 

L .388 , .0893 , 

6 = .083 deg, q = .029 deg sec , u 


.0197, .104 
.261 , -22.5 

-1 ~ 

= .026 msec , x 



.072 m . 


27 


(4.26) 



It is not always possible to find a satisfactory steady-state filter 
by these techniques, even when the system is completely observable. If 
there are any stable (or neutrally stable) modes of the open-loop system 
that are "uncontrollable" by the disturbances, then the steady-state 
variance and covariances of these modes will be zero, i.e. , the optimum 
time- varying filter will have estimated them exactly by the time it reaches 
the statistical steady-state. As a result, the steady-state gains for this 
mode will be zero; this may produce an unsatisfactory filter, particularly 
if the undisturbable mode is neutrally stable or only marginally stable 
(this will show up as a zero or very small eigenvalue for the steady-state 
filter). A straightforward remedy for this situation is to divide the 
system into disturbable and undisturbable modes, and then (a) design a KB 
filter for the disturbable modes neglecting the coupling to the undisturbr 
able modes (i.e., treat them as known forcing functions); then (b) design 
an "observer" for the undisturbable modes, neglecting the coupling to the 
disturbable modes, with eigenvalues very much less than the eigenvalues of 
the filter in (a), but large enough so that errors in these modes will die 
out rapidly enough to be satisfactory [see, e.g. , Bryson and Luenberger 
(1970) for observer design]. 

The following example illustrates a neutrally stable, undisturbable 
mode. An inverted pendulum is to be stabilized by torquing a wheel mounted 
at the tip of the pendulum as shown in the sketch below, A random wind 
blows on the pendulum rod but produces no torque on the wheel, so the 
system model is reasonably well described by 

"Undisturbable" would be a better word than "uncontrollable" in this 
context. 
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b = q 

q=n0+u+w 
r = -au 



where r = angular velocity of wheel, u is proportional to the torque 
applied to the wheel, and w is proportional to the wind velocity. Clearly, 
the wheel is undlsturbable by the wind but r is observable if r or r - q 
is measured; in principle, then, a nonstationary filter could bring r - r 
to zero exactly after a long time even with a very noisy measurement. The 
steady-state of such a filter would have zero gain on the measurements in 
the r channel, i.e., no attention would be paid to new measurement data. 


V. MEAN SQUARE BEHAVIOR OF SYSTEMS WITH FILTER/ STATE -FEEDBACK CONTROLLERS 
As discussed in the introduction, the control system that minimizes 
the quadratic criterion (2.3) on the average , in the presence of random 
disturbances, is constructed by using the state-feedback control gains 
of (2,5) on the estimated state from the filter (4.2), i.e,, 

u = -Cx (5.1) 

This combination of filter and state-feedback gains (see Fig. 1.2) consti- 
tutes a "compensator " lathe language of classical controller design. 

The system, controlled by this compensator, with random forcing 
functions, w(t), random errors in the sensor signals, v(t), and random 
initial conditions, x(t^), is itself a gauss-markov process. Hence, its mean 
and mean-square behavior can be predicted using (3.9) and (3.10), by considering 
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the state vector to be 


and the disturbance vector to be 


X is thus treated as the "state" of the compensator. However, this 


as 


prediction is actually simpler than just indicated if we use 
the state vector, where X := x - x = filter estimate error, because X 
does not depend on x : 

X 


(F - KH)X + Kv - rw ; X(t ) = -x(t ) 

o o 


X = (F - GC)x - KHX + Kv ; x(t ) = 0 

o 


where we have used 

z = Hx + v, 
u = -Cx 

fp 

Since E f [X(t )][x(t )] }=0, it follows from (5.2) that 
o o 


(5.2) 


E { CX(t)] Cx(t)]^) = 0 for t > t 


(5.3) 


(5.4) 


Therefore 


E{ [x(t)] [x(t)^} s E ( (x - x)(x - x) = E(xx^) + E(XX"^), 




or 


where 


X = X + P 


X = E(xx^), X - E(xx'^), P = E(Xx^) 


(5.5) 


Now P is given by (4.3) and it follows directly from (3.8), (5.2), 
and (5.4) that 

k = (F- CC)9. + S(F-(X:)^ +PH^R"^HP ; S(t ) =0 , (5.6) 

o 

T -1 

which is a linear matrix equation, since PH R HP is a known forcing 
function [i.e. , it can be pre-computed from (4.3)]. Fran (5.3) and (5.5) 
it follows that 

E[u(t)uT(t)] = CXC'T (5.7) 

For the statistically stationary case, P = 0, $ = 0, and (5.6) 

becomes a set of - 2 n(n + l) simultaneous linear equations. For even 
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moderate values of n, solution of that many linear equations is not 
a trivial task. 

For the helicopter example used in Sections II and IV, with the 
data of Tables 1 and 2, the predicted RMS state and control deviations 
in the gusty wind (RMS 7 msec”^, correlation time 5 sec.) using the 
state-feedback filter compensator are given in Table 4. 

To show the importance of both sensors in the example, we designed 
a controller using only the position sensor (it is not possible, of 
course, to use only the pitch angle sensor, since then the position 
error can not be estimated). For the same position sensor accuracy and 
the same gusty wind, the predicted RMS deviations were 9 = 24.4 deg., 

X = 4.76 m, 5 = 10.4 deg. which would be quite unacceptable. This 
shows that the pitch angle sensor is essential. 

The overall compensator transfer function , from the measurements, 
z(s), to the controls, u(s), for the stationary case is given by 

u(s) = -C[sl - F +GC +KH]-1 kz(s) (5.8) 

Some of the poles and zeroes of this transfer function may have positive 
real parts, i.e., the compensator does not necessarily have "minimum 
phase". (Most classical compensator designs have minimum phase. ) Thus, 
the compensator (5.8) will, in general, require active logic components, 
i.e., it cannot be built just from passive components like resistors and 
capacitors. However, the eigenvalues (poles) of the controlled system 
are simply the eigenvalues of the state feedback controller, F-GC, plus 
the eigenvalues of the Filter, F-KH, all of which have negative real 
parts. * 

♦ 

Thus, (5.8) may obscure the stability of the controlled system. 
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VI . SUMMARY 


The current theory available for designing controllers for 
dynamic systems subjected to random disturbances may be summarized 
as follows; 

(a) Model the system by a set of coupled first-order ordinary 
differential equations (i.e., put them in state vector form). 

(b) Linearize the system about a nominal path or operating 
condition. 

(c) Design a linear state-feedback controller using weighting 
factors in a quadratic performance index as described in 
Section II. 

(d) Model the random disturbances as Gauss-Markov processes, 
using shaping filters if necessary. 

(e) Design a Kalman-Bucy filter to estimate the state vector 
from the available measurements as described in Section IV. 

(f) Cascade the filter and the state-feedback controller, i.e., 
use state-feedback on the estimated state vector. This is 
the optimal "compensator" for the system. 

(g) Predict the RMS (root mean square) response of the controlled 
system in the presence of the random disturbances using the 
methods of Section V. If the response is unsatisfactory, 
modify the compensator by changing weighting factors in (c) 
and/or use additional or more accurate measurements. 
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FIG. 1. 


Block Diagram; Open- Loop (Programmed) Control Has No Feedback Loop 
Closed-Loop Control Has Feedback Loop 
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FIG. 2 


Block Diagram; Dynamic Compensator Formed 
from Filter and State-Feedback Gains 
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Nomenclature and Reference Axes 
for Helicopter Hover Example 
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FIG. 4. 


Root Locus for Position Regulator, Helicopter in Hover, 
with as Parameter 
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FIG. 8. 


Pendulum in Random Wing - Estimate Error Variance and Covariance 
Histories Using Measurement of q with Additive White Noise 
Having Spectral Density R 






FIG. 9. 


Root Locus for (x,6) Measurement Filter, Helicopter in Hover, 
with Q/R0 as Parameter, (R / ) = 5.6 X 10"5 m-2 

X V 
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TABLE 1 


REGULATOR DESIGN FOR A TYPICAL TRANSPORT HELICOPTER 


Coefficients: 

II 

.415 sec“^ , 

o = .0198 sec”^ , 

A 

, n = 6. 

27 sec 


a = 

. 0111 m“^ sec , 

a = 1.43 rnsec"^ , 

. g = 9. 

80 msec 


Open-Loop Eigenvalues and Eigenvectors* 


s, sec”^ 

0 

-.68 

+. 12 ±.38j 

0, deg. 

0 

2.4 

.79 ±.56j 

X, meters 

1.0 

1.0 

1.0 


Closed-Loop Eigenvalues and Eigenvectors* for (A^/B) = .002 m 


-1 

s , sec 

-1.24 ±.56j 

-.42 ±1. 14j 

6, deg. 

8.0 ^4.9j 

-5.4 T8. Ij 

X, meters 

1.0 

1.0 

--- - 


Feedback Gains for (A /B) = .002 m 

X 


C^ = .864 


C = .0869 m“l sec. 

C = .0447 m“l 

9 


u 

X 


where 6 = -Cq-Cu-Cx 

9 q u X 


The q and u components of the eigenvectors are given by q = s0 , 
u = sx. 
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TABLE 2 


FILTER DESIGN FOR A TYPICAL TRANSPORT HELICOPTER 
[RMS wind 7 msec~^ with correl. time 5 sec, Rx=.272 m^sec, Rg= 1 . 52 x lO'^sec] 

Eigenvalues of the Filter Estimate Errors (sec~^) 


-2.34 


-0.19 ± 0.19J 


-1.17 ± 2.00j 


Filter Gains (Units in meters, sec, radians) 


n 


q 

u 

X 

w 

K 

nx 



.000327 

.000251 

.0735 

.381 

-.128 

*^n0 

4.05 

8.21 

1 

21.4 

5.94 

-954 


RMS Estimate Errors 


0 = 0.45 deg 

q = 1.15 dog sec"^ 


QQQIHQI 

w = 3.78 msec“^ 
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J3> 


TABLE 3 


DYNAMIC COMPENSATOR FOR A TYPICAL TRANSPORT HELICOPTER 
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y\ 
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V- 

C u 
u 

- C X 
X 

- 

C 

w 

w = commanded rotor 

tilt , 



z = measurement of x, 

X 

Zg = measurement of 6, and 

C , K , K ^ given in Tables 1 and 2; = 5 sec, C = -.00195 m“l sec 

n nx n0 c ’ w 

TABLE 4 

PREDICTED RMS STATE AND CONTROL DEVIATIONS OF TYPICAL TRANSPORT HELICOPTER 
IN A GUSTY WIND USING CONTROLLER OF TABLES 1 AND 2 


a = 1.95 deg 

q = 2.92 deg sec“^ 

u = 0.25 msec"^ 

X =0.51 m 


5 = 1.33 deg 




